home
***
CD-ROM
|
disk
|
FTP
|
other
***
search
/
Collection of Tools & Utilities
/
Collection of Tools and Utilities.iso
/
fortran
/
linpkdrv.zip
/
SG.FOR
< prev
next >
Wrap
Text File
|
1984-01-05
|
18KB
|
616 lines
C MAIN PROGRAM
INTEGER LUNIT
C ALLOW 5000 UNDERFLOWS.
C CALL TRAPS(0,0,5001,0,0)
C
C OUTPUT UNIT NUMBER
C
LUNIT = 6
C
CALL SGETS(LUNIT)
C
STOP
END
SUBROUTINE SGETS(LUNIT)
C LUNIT IS THE OUTPUT UNIT NUMBER
C
C TESTS
C SGECO,SGEFA,SGESL,SGEDI,SGBCO,SGBFA,SGBSL,SGBDI
C
C LINPACK. THIS VERSION DATED 08/14/78 .
C CLEVE MOLER, UNIVERSITY OF NEW MEXICO, ARGONNE NATIONAL LAB.
C
C SUBROUTINES AND FUNCTIONS
C
C LINPACK SGECO,SGESL,SGEDI,SGBCO,SGBSL,SGBDI
C EXTERNAL SGEXX,SMACH
C BLAS SAXPY,SDOT,SSCAL,SASUM
C FORTRAN ABS,AMAX1,FLOAT,MAX0,MIN0
C
C INTERNAL VARIABLES
C
REAL A(15,15),AB(43,15),AINV(15,15),ASAVE(15,15)
REAL B(15),BT(15),SDOT,DET(2),DETB(2)
REAL X(15),XB(15),XEXACT(15),XT(15),XTB(15),T,Z(15)
REAL AINORM,ANORM,SMACH,COND,COND1,EN,ENORM,EPS
REAL ETNORM,FNI,FNORM,ONEPX,RCOND,RCONDB,RNORM
REAL RTNORM,Q(8),QS(8),SASUM,XNORM,XTNORM
INTEGER I,IPVT(15),IPVTB(15),IQ(8),I1,I2,J
INTEGER K,KASE,KB,KBFAIL,KOUNT,KP1,KSING,KSUSP(8)
INTEGER L,LDA,LDAB,LUNIT,M,ML,MU,N,NM1,NPRINT
LOGICAL KBF
C
LDA = 15
LDAB = 43
C
C WRITE MATRIX AND SOLUTIONS IF N .LE. NPRINT
C
NPRINT = 3
C
WRITE (LUNIT,460)
WRITE (LUNIT,880)
C
DO 10 I = 1, 8
KSUSP(I) = 0
10 CONTINUE
KSING = 0
KBFAIL = 0
C
C SET EPS TO ROUNDING UNIT
C
EPS = SMACH(1)
WRITE (LUNIT,470) EPS
WRITE (LUNIT,450)
C
C START MAIN LOOP
C
KASE = 1
20 CONTINUE
C
C GENERATE TEST MATRIX
C
CALL SGEXX(A,LDA,N,KASE,LUNIT)
C
C N = 0 SIGNALS NO MORE TEST MATRICES
C
C ...EXIT
IF (N .LE. 0) GO TO 440
ANORM = 0.0E0
DO 30 J = 1, N
ANORM = AMAX1(ANORM,SASUM(N,A(1,J),1))
30 CONTINUE
WRITE (LUNIT,650) ANORM
C
IF (N .GT. NPRINT) GO TO 50
WRITE (LUNIT,450)
DO 40 I = 1, N
WRITE (LUNIT,700) (A(I,J), J = 1, N)
40 CONTINUE
WRITE (LUNIT,450)
50 CONTINUE
C
C GENERATE EXACT SOLUTION
C
XEXACT(1) = 1.0E0
IF (N .GE. 2) XEXACT(2) = 0.0E0
IF (N .LE. 2) GO TO 70
DO 60 I = 3, N
XEXACT(I) = -XEXACT(I-2)
60 CONTINUE
70 CONTINUE
C
C SAVE MATRIX AND GENERATE R.H.S.
C
DO 90 I = 1, N
B(I) = 0.0E0
BT(I) = 0.0E0
DO 80 J = 1, N
ASAVE(I,J) = A(I,J)
B(I) = B(I) + A(I,J)*XEXACT(J)
BT(I) = BT(I) + A(J,I)*XEXACT(J)
80 CONTINUE
X(I) = B(I)
XT(I) = BT(I)
XB(I) = X(I)
XTB(I) = XT(I)
90 CONTINUE
C
C FACTOR AND ESTIMATE CONDITION
C
CALL SGECO(A,LDA,N,IPVT,RCOND,Z)
C
C OUTPUT NULL VECTOR IF N .LE. NPRINT
C
IF (N .GT. NPRINT) GO TO 110
WRITE (LUNIT,720)
DO 100 I = 1, N
WRITE (LUNIT,730) Z(I)
100 CONTINUE
WRITE (LUNIT,450)
110 CONTINUE
C
C FACTOR BAND FORM AND COMPARE
C
KBF = .FALSE.
ML = 0
MU = 0
DO 140 J = 1, N
DO 130 I = 1, N
IF (ASAVE(I,J) .EQ. 0.0E0) GO TO 120
IF (I .LT. J) MU = MAX0(MU,J-I)
IF (I .GT. J) ML = MAX0(ML,I-J)
120 CONTINUE
130 CONTINUE
140 CONTINUE
WRITE (LUNIT,790) ML,MU
IF (2*ML + MU + 1 .LE. LDAB) GO TO 150
WRITE (LUNIT,680)
GO TO 430
150 CONTINUE
M = ML + MU + 1
DO 170 J = 1, N
I1 = MAX0(1,J-MU)
I2 = MIN0(N,J+ML)
DO 160 I = I1, I2
K = I - J + M
AB(K,J) = ASAVE(I,J)
160 CONTINUE
170 CONTINUE
C
CALL SGBCO(AB,LDAB,N,ML,MU,IPVTB,RCONDB,Z)
C
IF (RCONDB .EQ. RCOND) GO TO 180
WRITE (LUNIT,780)
WRITE (LUNIT,820) RCOND,RCONDB
KBF = .TRUE.
180 CONTINUE
KOUNT = 0
DO 190 J = 1, N
IF (AB(M,J) .NE. A(J,J)) KOUNT = KOUNT + 1
IF (IPVTB(J) .NE. IPVT(J)) KOUNT = KOUNT + 1
190 CONTINUE
IF (KOUNT .EQ. 0) GO TO 200
WRITE (LUNIT,780)
WRITE (LUNIT,830) KOUNT
KBF = .TRUE.
200 CONTINUE
C
C TEST FOR SINGULARITY
C
IF (RCOND .GT. 0.0E0) GO TO 210
WRITE (LUNIT,710) RCOND
WRITE (LUNIT,480)
KSING = KSING + 1
GO TO 420
210 CONTINUE
COND = 1.0E0/RCOND
WRITE (LUNIT,500) COND
ONEPX = 1.0E0 + RCOND
IF (ONEPX .EQ. 1.0E0) WRITE (LUNIT,490)
C
C COMPUTE INVERSE, DETERMINANT AND COND1 = TRUE CONDITION
C
DO 230 J = 1, N
DO 220 I = 1, N
AINV(I,J) = A(I,J)
220 CONTINUE
230 CONTINUE
CALL SGEDI(AINV,LDA,N,IPVT,DET,Z,11)
AINORM = 0.0E0
DO 240 J = 1, N
AINORM = AMAX1(AINORM,SASUM(N,AINV(1,J),1))
240 CONTINUE
COND1 = ANORM*AINORM
WRITE (LUNIT,510) COND1
WRITE (LUNIT,750) DET(1)
WRITE (LUNIT,760) DET(2)
C
C SOLVE A*X = B AND TRANS(A)*XT = BT
C
CALL SGESL(A,LDA,N,IPVT,X,0)
CALL SGESL(A,LDA,N,IPVT,XT,1)
C
IF (N .GT. NPRINT) GO TO 270
WRITE (LUNIT,520)
DO 250 I = 1, N
WRITE (LUNIT,740) X(I)
250 CONTINUE
WRITE (LUNIT,530)
DO 260 I = 1, N
WRITE (LUNIT,740) XT(I)
260 CONTINUE
WRITE (LUNIT,450)
270 CONTINUE
C
C MORE BAND COMPARE
C
CALL SGBSL(AB,LDAB,N,ML,MU,IPVTB,XB,0)
CALL SGBSL(AB,LDAB,N,ML,MU,IPVTB,XTB,1)
KOUNT = 0
DO 280 I = 1, N
IF (XB(I) .NE. X(I)) KOUNT = KOUNT + 1
IF (XTB(I) .NE. XT(I)) KOUNT = KOUNT + 1
280 CONTINUE
IF (KOUNT .EQ. 0) GO TO 290
WRITE (LUNIT,780)
WRITE (LUNIT,840) KOUNT
KBF = .TRUE.
290 CONTINUE
CALL SGBDI(AB,LDAB,N,ML,MU,IPVTB,DETB)
IF (DETB(1) .EQ. DET(1) .AND. DETB(2) .EQ. DET(2))
* GO TO 300
WRITE (LUNIT,780)
WRITE (LUNIT,850) DETB
KBF = .TRUE.
300 CONTINUE
C
C RECONSTRUCT A FROM TRIANGULAR FACTORS , L AND U
C
NM1 = N - 1
IF (NM1 .LT. 1) GO TO 330
DO 320 KB = 1, NM1
K = N - KB
KP1 = K + 1
L = IPVT(K)
DO 310 J = KP1, N
T = -A(K,J)
CALL SAXPY(N-K,T,A(K+1,K),1,A(K+1,J),1)
T = A(L,J)
A(L,J) = A(K,J)
A(K,J) = T
310 CONTINUE
T = -A(K,K)
CALL SSCAL(N-K,T,A(K+1,K),1)
T = A(L,K)
A(L,K) = A(K,K)
A(K,K) = T
320 CONTINUE
330 CONTINUE
C
C COMPUTE ERRORS AND RESIDUALS
C E = X - XEXACT
C ET = XT - XEXACT
C R = B - A*X
C RT = BT - A*XT
C F = A - L*U
C AI = A*INV(A) - I
C
XNORM = SASUM(N,X,1)
XTNORM = SASUM(N,XT,1)
ENORM = 0.0E0
ETNORM = 0.0E0
FNO